home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / linnew.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  26.2 KB  |  899 lines
  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  8. ;;;     (c) Copyright 1980 Massachusetts Institute of Technology         ;;;
  9. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10.  
  11. (in-package "MAXIMA")
  12. (macsyma-module linnew)
  13.  
  14. ;; This is a matrix package which uses minors, basically.
  15. ;; TMLINSOLVE(LIST-OF-EQUAIONS,LIST-OF-VARIABLES,LIST-OF-VARIABLES-TO-BE-OBTAINED)
  16. ;; solves the linear equation. LIST-OF-VARIABLES-TO-BE-OBTAINED can be omitted,
  17. ;; in which case all variables are obtained. TMNEWDET(MATRIX,DIMENSION)
  18. ;; computes the determinant.  DIMENSION can be omitted.  The default is
  19. ;; DIMENSION=(declared dimension of MATRIX). TMINVERSE(MATRIX) computes the
  20. ;; inverse of matrix.
  21.  
  22. ;; The program uses hash arrays to remember the minors if N > threshold.  If
  23. ;; $WISE is set to T, the program knocks out unnecessary elements.  But also it
  24. ;; kills necessary ones in the case of zero elements! The $WISE flag should
  25. ;; not be set to T for inverse.  The default of $WISE is NIL.
  26.  
  27. ;; There seem to have been a number of bugs in this code.  I changed
  28. ;; the array referencing to value cell, and some of the stuff about
  29. ;; cre form.  It now seems tminverse  and tmlinsolve, now seem to work. --wfs.
  30.  
  31. ;;these are arrays
  32. (DECLARE-TOP(special  *TMARRAYS*  *A2*  *B*  *AA* 
  33.              *ROW*  *COL*  *ROWINV*  *COLINV*  *INDX* ))
  34.  
  35. (DECLARE-TOP(SPECIAL N NX IX)) 
  36.  
  37. (DECLARE-TOP(SPECIAL $LINENUM $DISPFLAG $LINECHAR $WISE $FOOL)) 
  38.  
  39. (defvar *tmarrays* nil)
  40.  
  41. ;; If N < threshold declared array is used, otherwise hashed array.
  42.  
  43.  
  44. (DEFMACRO THRESHOLD () 10.)
  45.  
  46. (DEFUN TMINITIALFLAG NIL 
  47.        (COND ((NOT (BOUNDP '$WISE)) (SETQ $WISE NIL)))
  48.        (COND ((NOT (BOUNDP '$FOOL)) (SETQ $FOOL NIL))))
  49.  
  50. ;; TMDET returns the determinant of N*N matrix A2 which is in an globally
  51. ;; declared array A2.
  52.  
  53. (DEFUN TMDET (A4 N) 
  54.        (PROG (INDEX RESULT IX) 
  55.          (TMINITIALFLAG)
  56.          (TMHEADING)
  57.          (SETQ IX 0. NX 0.)
  58.          (DO ((I 1. (f1+ I)))
  59.          ((> I N))
  60.          (SETQ INDEX (CONS I INDEX)))
  61.          (SETQ INDEX ;(REVERSE INDEX)
  62.                  (nreverse index))
  63.          (SETQ RESULT (TMINOR A4 N 1. INDEX 0.))
  64.          (RETURN RESULT)))
  65.  
  66. ;; TMLIN SOLVES M SETS OF LINEAR EQUATIONS WHITH N UNKNOWN VARIABLES. IT SOLVES
  67. ;; ONLY FOR THE FIRST NX UNKNOWNS OUT OF N. THE EQUATIONS ARE EXPRESSED IN
  68. ;; MATRIX FORM WHICH IS IN N*(N+M) ARRAY A2. AS USUAL , THE LEFT HAND SIDE N*N
  69. ;; OF A2 REPRESENTS THE COEFFICIENT MATRIX, AND NEXT N*M OF A2 IS THE RIGHT
  70. ;; HAND SIDE OF THE M SETS OF EQUATIONS.  SUPPOSE N=3, M=2, AND THE UNKKNOWNS
  71. ;; ARE (X1 Y1 Z1) FOR THE FIRST SET AND (X2 Y2 Z2) FOR THE SECOND. THEN THE
  72. ;; RESULT OF TMLIN IS ((DET) (U1 U2) (V1 V2) (W1 W2)) WHERE DET IS THE
  73. ;; DETERMINANT OF THE COEFFICIENT MATRIX AND X1=U1/DET, X2=U2/DET, Y1=V1/DET,
  74. ;; Y2=V2/DET ETC.
  75.  
  76. (DEFUN TMLIN (A4 N M NX) 
  77.        (PROG (INDEX R) 
  78.          (TMDEFARRAY N)
  79.          (TMINITIALFLAG)
  80.          (TMHEADING)
  81.          (DO ((I 1. (f1+ I)))
  82.          ((> I N))
  83.          (SETQ INDEX (CONS I INDEX)))
  84.          (SETQ INDEX (REVERSE INDEX))
  85.          (SETQ R
  86.            (DO ((IX 0. (f1+ IX)) (RESULT))
  87.                ((> IX NX) (REVERSE RESULT))
  88.                (SETQ RESULT
  89.                  (CONS (DO ((I 1. (f1+ I)) (RES))
  90.                        ((> I
  91.                        (COND ((= IX 0.) 1.)
  92.                          (T M)))
  93.                     (REVERSE RES))
  94.                        (COND ((NOT $WISE)
  95.                           (TMKILLARRAY IX)))
  96.                        (SETQ RES
  97.                          (CONS (TMINOR A4
  98.                                N
  99.                                1.
  100.                                INDEX
  101.                                I)
  102.                            RES)))
  103.                    RESULT))
  104.                (COND ((AND (= IX 0.)
  105.                    (EQUAL (CAR RESULT)
  106.                       '(0. . 1.)))
  107.                   (merror "COEFFICIENT MATRIX IS SINGULAR")))))
  108.          (TMREARRAY N)
  109.          (RETURN R)))
  110.  
  111. ;; TMINOR ACTUALLY COMPUTES THE MINOR DETERMINANT OF A SUBMATRIX OF A2, WHICH
  112. ;; IS CONSTRUCTED BY EXTRACTING ROWS (K,K+1,K+2,...,N) AND COLUMNS SPECIFIED BY
  113. ;; INDEX. N IS THE DIMENSION OF THE ORIGINAL MATRIX A2.  WHEN TMINOR IS USED
  114. ;; FOR LINEAR EQUATION PROGRAM, JRIGHT SPECIFIES A COLUMN OF THE CONSTANT
  115. ;; MATRIX WHICH IS PLUGED INTO AN IX-TH COLUMN OF THE COEFFICIENT MATRIX FOR
  116. ;; ABTAINING IX-TH UNKNOWN. IN OTHER WORDS, JRIGHT SPECIFIES JRIGHT-TH
  117. ;; EQUATION.
  118.  
  119.  
  120. (DEFUN TMINOR (A4 N K INDEX JRIGHT) 
  121.        (PROG (SUBINDX L RESULT NAME AORB)
  122.          (setq a4 (get-array-pointer a4))
  123.          (COND
  124.           ((= K N)
  125.            (SETQ RESULT
  126.              (COND ((= K IX) (AREF A4 (CAR INDEX) (f+ JRIGHT N)))
  127.                (T (AREF A4 (CAR INDEX) K)))))
  128.           (T
  129.            (DO
  130.         ((J 1. (f1+ J)) (SUM '(0. . 1.)))
  131.         ((> J (f1+ (f- N K))) (SETQ RESULT SUM))
  132.         (SETQ L (EXTRACT INDEX J))
  133.         (SETQ SUBINDX (CADR L))
  134.         (SETQ L (CAR L))
  135.         (SETQ AORB (COND ((= K IX) (AREF A4 L (f+ JRIGHT N)))
  136.                  (T (AREF A4 L K))))
  137.         (COND
  138.          ((NOT (EQUAL AORB '(0. . 1.)))
  139.           (SETQ NAME (TMACCESS SUBINDX))
  140.           (SETQ 
  141.            SUM
  142.            (funcall (COND ((ODDP J) 'RATPLUS)
  143.                   (T 'RATDIFFERENCE))
  144.             SUM
  145.             (RATTIMES
  146.              AORB
  147.              (COND ($FOOL (TMINOR A4 N (f1+ K) SUBINDX JRIGHT))
  148.                (T (COND ((NOT (NULL (TMEVAL NAME)))
  149.                      (TMEVAL NAME))
  150.                     ((TMNOMOREUSE J L K)
  151.                      (TMSTORE NAME NIL)
  152.                      (TMINOR A4
  153.                          N
  154.                          (f1+ K)
  155.                          SUBINDX
  156.                          JRIGHT))
  157.                     (T (TMSTORE NAME
  158.                         (TMINOR A4
  159.                             N
  160.                             (f1+ K)
  161.                             SUBINDX
  162.                             JRIGHT))))))
  163.              T)))))
  164.         (COND ($WISE (COND ((TMNOMOREUSE J L K)
  165.                     (TMKILL SUBINDX K))))))))
  166.          (RETURN RESULT))) 
  167.  
  168. (DEFUN EXTRACT (INDEX J) 
  169.        (DO ((IND INDEX (CDR IND)) (COUNT 1. (f1+ COUNT)) (SUBINDX))
  170.        ((NULL IND))
  171.        (COND ((= COUNT J)
  172.           (RETURN (LIST (CAR IND) (NCONC SUBINDX (CDR IND)))))
  173.          (T (SETQ SUBINDX (NCONC SUBINDX (LIST (CAR IND)))))))) 
  174.  
  175. (DECLARE-TOP(SPECIAL VLIST VARLIST GENVAR)) 
  176.  
  177.  
  178. (DEFUN TMRATCONV (BBB N M) 
  179.        (PROG (CCC)
  180.          (declare (special ccc))    ;Tell me this worked in Maclisp.  --gsb
  181.          ;Actually, i suspect it didn't, at least ever since
  182.          ; (sstatus punt).
  183.          (SET 'CCC BBB)
  184.          (DO ((K 1. (f1+ K)))
  185.          ((> K N))
  186.          (DO ((J 1. (f1+ J)))
  187.              ((> J M))
  188.              (NEWVAR1 (STORE (aref *A2* K J)
  189.                      (maref ccc k j)
  190. ;;                     (nth j (nth k *a2*))
  191. ;;                     (MEVAL (LIST (LIST 'CCC 'array) K J))  ;;just the
  192.                      ))))
  193.          
  194.          (NEWVAR (CONS '(MTIMES) VLIST))
  195.          (DO ((K 1. (f1+ K)))
  196.          ((> K N))
  197.          (DO ((J 1. (f1+ J)))
  198.              ((> J M))
  199.              (STORE (aref *A2* K J)
  200.                 (CDR (RATREP* (aref *A2* K J)))))))) 
  201.  
  202. (DEFMFUN $TMNEWDET N 
  203.        (PROG (*AA* R VLIST) 
  204.          (COND ((= N 2.)
  205.             (COND ((NOT (INTEGERP (SETQ N (ARG 2.))))
  206.                (merror  "WRONG ARG")))
  207.             (SETQ *AA* (ARG 1.)))
  208.            ((AND (= N 1.) ($MATRIXP (SETQ *AA* (ARG 1.))))
  209.             (SETQ N (LENGTH (CDR (ARG 1.)))))
  210.            (T (merror "WRONG ARG")))
  211.          (setq  *A2* (*ARRAY nil 'T (f1+ N) (f1+ N)))
  212.          (TMDEFARRAY N)
  213.          (TMRATCONV *AA* N N)
  214.          (SETQ R (CONS (LIST 'MRAT
  215.                  'SIMP
  216.                  VARLIST
  217.                   GENVAR)
  218.                (TMDET '*A2* N)))
  219.          (*TMREARRAY '*A2*)
  220.          (TMREARRAY N)
  221.          (RETURN R))) 
  222.  
  223. (DEFMFUN $TMLINSOLVE NARG (TMLINSOLVE (LISTIFY NARG))) 
  224.  
  225. (DEFUN TMLINSOLVE (ARGLIST) 
  226.        (PROG (EQUATIONS VARS OUTVARS RESULT *AA*) 
  227.          (SETQ EQUATIONS (CDAR ARGLIST) 
  228.            VARS (CDADR ARGLIST) 
  229.            OUTVARS (COND ((NULL (CDDR ARGLIST)) VARS)
  230.                  (T (CDADDR ARGLIST))) 
  231.            ARGLIST NIL)
  232.          (SETQ VARS (TMERGE VARS OUTVARS))
  233.          (SETQ NX (LENGTH OUTVARS))
  234.          (SETQ N (LENGTH VARS))
  235.          (COND ((NOT (= N (LENGTH EQUATIONS)))
  236.             (RETURN (PRINT 'TOO-FEW-OR-MUCH-EQUATIONS))))
  237.          (SETQ 
  238.           *AA*
  239.           (CONS
  240.            '($MATRIX SIMP)
  241.            (MAPCAR 
  242.         #'(LAMBDA (EXP) 
  243.           (APPEND
  244.            '((MLIST))
  245.            (MAPCAR #'(LAMBDA (V) 
  246.                     (PROG (R) 
  247.                       (SETQ EXP
  248.                         ($BOTHCOEF EXP V)
  249.                         R
  250.                         (CADR EXP)
  251.                         EXP
  252.                         (MEVAL (CADDR EXP)))
  253.                       (RETURN R)))
  254.                VARS)
  255.            (LIST (LIST '(MMINUS) EXP))))
  256.         (MAPCAR #'(LAMBDA (E) (MEVAL (LIST '(MPLUS)
  257.                           ($LHS E)
  258.                           (LIST '(MMINUS)
  259.                             ($RHS E)))))
  260.             EQUATIONS))))
  261.          (SETQ RESULT (CDR ($TMLIN *AA* N 1. NX)))
  262.          (RETURN
  263.           (DO
  264.            ((VARS (CONS NIL OUTVARS) (CDR VARS))
  265.         (LABELS)
  266.         (DLABEL)
  267.         (NAME))
  268.            ((NULL VARS)
  269.         (CONS '(MLIST) (CDR (REVERSE LABELS))))
  270.            (SETQ NAME (MAKELABEL $LINECHAR))
  271.            (SETQ $LINENUM (f1+ $LINENUM))
  272.            (SET NAME
  273.             (COND ((NULL (CAR VARS))
  274.                (SETQ DLABEL NAME)
  275.                (CADAR RESULT))
  276.               (T (LIST '(MEQUAL)
  277.                    (CAR VARS)
  278.                    (LIST '(MTIMES SIMP)
  279.                      (CADAR RESULT)
  280.                      (LIST '(MEXPT SIMP)
  281.                            DLABEL
  282.                            -1.))))))
  283.            (SETQ LABELS (CONS NAME LABELS))
  284.            (SETQ RESULT (CDR RESULT))
  285.            (COND
  286.         ($DISPFLAG (MTELL-OPEN "~M" (NCONC (NCONS '(MLABLE))
  287.                       (NCONS NAME)
  288.                       (NCONS (EVAL NAME)))))))))) 
  289.  
  290. (DEFUN TMERGE (VARS OUTVARS) 
  291.        (APPEND OUTVARS
  292.            (PROG (L) 
  293.              (MAPCAR #'(LAMBDA (V) 
  294.                       (COND ((zl-MEMBER V OUTVARS) NIL)
  295.                         (T (SETQ L (CONS V L)))))
  296.                  VARS)
  297.              (RETURN (REVERSE L))))) 
  298.  
  299. (DEFMFUN $TMLIN (*AA* N M NX) 
  300.        (PROG (R VLIST) 
  301.          (setq  *A2* (*ARRAY nil 'T (f1+ N) (f1+ (f+ M N))))
  302.          (show *a2*)
  303.          (TMRATCONV *AA* N (f+ M N))
  304.          (SETQ 
  305.           R
  306.           (CONS
  307.            '(MLIST)
  308.            (MAPCAR 
  309.         #'(LAMBDA (RES) 
  310.           (CONS '(MLIST)
  311.             (MAPCAR #'(LAMBDA (RESULT) 
  312.                      (CONS (LIST 'MRAT
  313.                              'SIMP
  314.                              VARLIST
  315.                               GENVAR)
  316.                            RESULT))
  317.                 RES)))
  318.         (TMLIN '*A2* N M NX))))
  319.          (*TMREARRAY '*A2*)
  320.          (show *a2*)
  321.          (RETURN R))) 
  322.  
  323. (DEFUN TMKILL (*INDX* K) 
  324.        (PROG (NAME SUBINDX J L) 
  325.          (COND ((NULL *INDX*) (RETURN NIL)))
  326.          (SETQ NAME (TMACCESS *INDX*))
  327.          (COND ((NOT (NULL (TMEVAL NAME))) (TMSTORE NAME NIL))
  328.            (T (DO ((IND *INDX* (CDR IND)) (COUNT 1. (f1+ COUNT)))
  329.               ((NULL IND))
  330.               (SETQ L (EXTRACT *INDX* COUNT) 
  331.                 J (CAR L) 
  332.                 SUBINDX (CADR L))
  333.               (COND ((= J COUNT)
  334.                  (TMKILL SUBINDX (f1+ K))))))))) 
  335.  
  336. (DEFUN TMNOMOREUSE (J L K) 
  337.        (COND ((AND (= J L) (OR (> K NX) (< K (f1+ IX)))) T) (T NIL))) 
  338.  
  339. (DEFUN TMDEFARRAY (N) 
  340.        (PROG (NAME) 
  341.          (COND
  342.           (;(GET '*TMARRAYS* 'array)
  343.            (setq *tmarrays* (get-array-pointer *tmarrays*))
  344.            (TMREARRAY (f1- (COND ((CADR (ARRAYDIMS *TMARRAYS*)))
  345.                     (T 1.))))))
  346.          (setq  *TMARRAYS* (*ARRAY nil 'T (f1+ N)))
  347.          (DO ((I 1. (f1+ I)))
  348.          ((> I N))
  349.          (SETQ NAME (COND ((= I 1.) (make-symbol "M"))
  350.                   (T (GENSYM))))
  351.          (COND ((< N (THRESHOLD))
  352.             ;(STORE (aref *TMARRAYS* I) NAME)
  353.              (set name (*ARRAY nil T (f1+ (TMCOMBI N I))))
  354.              (STORE (aref *TMARRAYS* I) (get-array-pointer NAME))
  355.              )
  356.                
  357.                (T (STORE (aref *TMARRAYS* I)
  358.                  (LIST NAME
  359.                        'SIMP
  360.                        'array)))))
  361.          (GENSYM "G")))
  362.  
  363. ;; TMREARRAY kills the TMARRAYS which holds pointers to minors. If (TMARRAYS I)
  364. ;; is an atom, it is declared array.  Otherwise it is hashed array.
  365.  
  366. (DEFUN TMREARRAY (N) 
  367.        (PROG NIL 
  368.          (DO ((I 1. (f1+ I)))
  369.          ((> I N))
  370.          (COND ((ATOM (aref *TMARRAYS* I)) (*TMREARRAY (aref *TMARRAYS* I)))
  371.                (T (TM$KILL (CAR (aref *TMARRAYS* I))))))
  372.          (*TMREARRAY '*TMARRAYS*))) 
  373.  
  374. (DEFUN TMACCESS (INDEX) 
  375.        (PROG (L) 
  376.          (COND ($FOOL (RETURN NIL)))
  377.          (SETQ L (LENGTH INDEX))
  378.          (RETURN
  379.           (COND ((< N (THRESHOLD))
  380.              (LIST 'aref (aref *TMARRAYS* L)
  381.                (DO ((I 1. (f1+ I))
  382.                 (X 0. (CAR Y))
  383.                 (Y INDEX (CDR Y))
  384.                 (SUM 0.))
  385.                    ((> I L) (f1+ SUM))
  386.                    (DO ((J (f1+ X) (f1+ J)))
  387.                    ((= J (CAR Y)))
  388.                    (SETQ SUM (f+ SUM
  389.                         (TMCOMBI (f- N J)
  390.                              (f- L I))))))))
  391.             (T (cons 'aref (CONS (aref *TMARRAYS* L) INDEX)))))) )
  392.  
  393. (DEFUN TMCOMBI (N I) 
  394.        (COND ((> (f- N I) I)
  395.           (// (TMFACTORIAL N (f- N I)) (TMFACTORIAL I 0.)))
  396.          (T (// (TMFACTORIAL N I) (TMFACTORIAL (f- N I) 0.))))) 
  397.  
  398. (DEFUN TMFACTORIAL (I J) 
  399.        (COND ((= I J) 1.) (T (f* I (TMFACTORIAL (f1- I) J))))) 
  400.  
  401. (DEFUN TMSTORE (NAME X) 
  402.        (COND ((< N (THRESHOLD))
  403.           (EVAL (LIST 'STORE NAME (LIST 'QUOTE X))))
  404.          (T (MSET NAME (LIST '(MQUOTE SIMP) X)) X)))
  405.  
  406. ;; TMKILLARRAY kills all (N-IX+1)*(N-IX+1) minors which are not necessary for
  407. ;; the computation of IX-TH variable in the linear equation.  Otherwise, they
  408. ;; will do harm.
  409.  
  410. (DEFUN TMKILLARRAY (IX) 
  411.        (DO ((I (f1+ (f- N IX)) (f1+ I)))
  412.        ((> I N))
  413.        (COND ((< N (THRESHOLD))
  414.           (FILLARRAY (aref *TMARRAYS* I) '(NIL)))
  415.          (T (TM$KILL (CAR (aref *TMARRAYS* I))))))) 
  416.  
  417. (DEFUN TMHEADING NIL NIL) 
  418.  
  419. (DEFUN TMEVAL (E) 
  420.        (PROG (RESULT) 
  421.          (RETURN (COND ((< N (THRESHOLD)) (EVAL E))
  422.                (T (SETQ RESULT (MEVAL E))
  423.                   (COND ((EQUAL RESULT E) NIL)
  424.                     (T (CADR RESULT)))))))) 
  425.  
  426. (DEFUN TM$KILL (E) (KILL1 E))
  427.  
  428. (DEFMFUN $TMINVERSE ( *AA*) 
  429.        (PROG (R VLIST N M NX) 
  430.          (SETQ N (LENGTH (CDR *AA*)) M N NX N)
  431.          (setq  *A2* (*ARRAY nil 'T (f1+ N) (f1+ (f+ M N))))
  432.          (TMRATCONV *AA* N N)
  433.          (DO ((I 1. (f1+ I)))
  434.          ((> I N))
  435.          (DO ((J 1. (f1+ J)))
  436.              ((> J M))
  437.              (STORE (aref *A2* I (f+ N J))
  438.                 (COND ((= I J) '(1. . 1.))
  439.                   (T '(0. . 1.))))))
  440.          (SETQ 
  441.           R
  442.           (MAPCAR 
  443.            #'(LAMBDA (RES) 
  444.          (CONS
  445.           '(MLIST)
  446.           (MAPCAR 
  447.            #'(LAMBDA (RESULT) 
  448.                 ($RATDISREP (CONS (LIST 'MRAT
  449.                             'SIMP
  450.                             VARLIST
  451.                              GENVAR)
  452.                           RESULT)))
  453.            RES)))
  454.            (TMLIN '*A2* N M NX)))
  455.          (SETQ R
  456.            (LIST '(MTIMES SIMP)
  457.              (LIST '(MEXPT SIMP) (CADAR R) -1.)
  458.              (CONS '($MATRIX SIMP) (CDR R))))
  459.          (*TMREARRAY '*A2*)
  460.          (RETURN R))) 
  461.  
  462. (DEFUN *TMREARRAY (X) (*REARRAY X)) 
  463.  
  464. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;                   
  465. ;;THIS IS A UTILITY PACKAGE FOR SPARSE
  466. ;;MATRIX INVERSION. A3 IS A N*N MATRIX.
  467. ;;IT RETURNS A LIST OF LISTS, SUCH AS
  468. ;;((I1 I2 ...) (J1 J2...) ...) WHERE (I1
  469. ;;I2 ..) SHOWS THE ROWS WHICH BELONGS TO
  470. ;;THE FIRST BLOCK, AND SO ON.  THE ROWS
  471. ;;SHOUD BE REORDERED IN THIS ORDER. THE
  472. ;;COLUMNS ARE NOT CHANGED. IT RETURNS NIL
  473. ;;IF A3 IS "OBVIOUSLY" SINGULAR.
  474.  
  475. ;; (DEFUN TMISOLATE (A3 N)
  476. ;;        (PROG (NODELIST)
  477. ;;          (SETQ A3 (GET A3 'ARRAY))
  478. ;;          (setq  B (*ARRAY nil 'T (f1+ N) (f1+ N)))
  479. ;;          (setq  ROW (*ARRAY nil 'T (f1+ N)))
  480. ;;          (setq  COL (*ARRAY nil 'T (f1+ N)))
  481. ;;          (DO ((I 1. (f1+ I)))
  482. ;;          ((> I N))
  483. ;;          (STORE (ROW I) I)
  484. ;;          (STORE (COL I) I))
  485. ;;          (DO ((I 1. (f1+ I)))
  486. ;;          ((> I N))
  487. ;;          (DO ((J 1. (f1+ J)))
  488. ;;              ((> J N))
  489. ;;              (STORE (B I J)
  490. ;;                 (NOT (EQUAL (AREF A3 I J)
  491. ;;                     '(0. . 1.))))))
  492. ;;          (COND ((NULL (TMPIVOT-ISOLATE 1.))
  493. ;;             (SETQ NODELIST NIL)
  494. ;;             (GO EXIT)))
  495. ;;          (DO ((I 1. (f1+ I)))
  496. ;;          ((> I N))
  497. ;;          (DO ((J 1. (f1+ J)))
  498. ;;              ((> J I))
  499. ;;              (STORE (B (ROW J) (COL I))
  500. ;;                 (OR (B (ROW I) (COL J))
  501. ;;                 (B (ROW J) (COL I))))
  502. ;;              (STORE (B (ROW I) (COL J)) (B (ROW J) (COL I))))
  503. ;;          (STORE (B (ROW I) (COL I)) T))
  504. ;;          (DO ((I 1. (f1+ I)))
  505. ;;          ((> I N))
  506. ;;          (COND ((EQ (B (ROW I) (COL I)) T)
  507. ;;             (SETQ NODELIST
  508. ;;                   (CONS (TMPULL-OVER I N) NODELIST)))))
  509. ;;          EXIT
  510. ;;          (*TMREARRAY 'B)
  511. ;;          (*TMREARRAY 'ROW)
  512. ;;          (*TMREARRAY 'COL)
  513. ;;          (RETURN (REVERSE NODELIST))))) 
  514.  
  515. ;; (DEFUN TMPULL-OVER (P N) 
  516. ;;        (PROG (Q) 
  517. ;;          (STORE (B (ROW P) (COL P)) NIL)
  518. ;;          (DO ((J 1. (f1+ J)))
  519. ;;          ((> J N) (SETQ Q NIL))
  520. ;;          (COND ((EQ (B (ROW P) (COL J)) T)
  521. ;;             (RETURN (SETQ Q J)))))
  522. ;;          (COND ((NULL Q) (RETURN (LIST (ROW P))))
  523. ;;            (T (DO ((J 1. (f1+ J)))
  524. ;;               ((> J N))
  525. ;;               (STORE (B (ROW Q) (COL J))
  526. ;;                  (OR (B (ROW Q) (COL J))
  527. ;;                      (B (ROW P) (COL J))))
  528. ;;               (STORE (B (ROW J) (COL Q))
  529. ;;                  (B (ROW Q) (COL J))))
  530. ;;               (TMCRIP P)
  531. ;;               (RETURN (CONS (ROW P) (TMPULL-OVER Q N))))))) 
  532.  
  533. ;; (DEFUN TMCRIP (P) 
  534. ;;        (DO ((I 1. (f1+ I)))
  535. ;;        ((> I N))
  536. ;;        (STORE (B (ROW P) (COL I)) NIL)
  537. ;;        (STORE (B (ROW I) (COL P)) NIL)))        
  538.  
  539. ;;TMPIVOT-ISOLATE CARRIES OUT PIVOTTING
  540. ;;SO THAT THE ALL DIAGONAL ELEMENTS ARE
  541. ;;NONZERO. THIS GARANTIES WE HAVE MAXIMUM
  542. ;;NUMBER OF BLOCKS ISOLATED.
  543.  
  544. (DEFUN TMPIVOT-ISOLATE (K) 
  545.        (COND ((> K N) T)
  546.          (T (DO ((I K (f1+ I)))
  547.             ((> I N) NIL)
  548.             (COND ((aref *B* (aref *ROW* I) (aref *COL* K))
  549.                (TMEXCHANGE '*ROW* K I)
  550.                (COND ((TMPIVOT-ISOLATE (f1+ K)) (RETURN T))
  551.                  (T (TMEXCHANGE '*ROW*
  552.                         K
  553.                         I))))))))) 
  554.  
  555. (DEFUN TMEXCHANGE (ROWCOL I J) 
  556.        (PROG (DUMMY) 
  557.          (SETQ ROWCOL (get-array-pointer rowcol))
  558.          (SETQ DUMMY (AREF ROWCOL I))
  559.          (STORE (AREF ROWCOL I) (AREF ROWCOL J))
  560.          (STORE (AREF ROWCOL J) DUMMY)))    
  561.  
  562.  
  563. ;; PROGRAM TO PREDICT ZERO ELEMENTS IN
  564. ;; THE SOLUTION OF INVERSE OR LINEAR
  565. ;; EQUATION. A IS THE COEFFICIENT MATRIX.
  566. ;; B IS THE RIGHT HAND SIDE MATRIX FOR
  567. ;; LINEAR EQUATIONS. A3 IS N*N AND B IS
  568. ;; M*M. X IS AN N*M MATRIX WHERE T -NIL
  569. ;; PATTERN SHOWING THE ZERO ELEMENTS IN
  570. ;; THE RESULT IS RETURND. T CORRESPONDS TO
  571. ;; NON-ZERO ELEMENT. IN THE CASE OF
  572. ;; INVERSE, YOU CAN PUT ANYTHING (SAY,NIL)
  573. ;; FOR B AND 0 FOR M.  NORMALLY IT RETURNS
  574. ;; T, BUT IN CASE OF SINGULAR MATRIX, IT
  575. ;; RETURNS NIL.
  576.  
  577. ;; (DEFUN TMPREDICT (A3 B X N M)
  578. ;;   (PROG (FLAGINV FLAG-NONSINGULAR)
  579. ;;     (SETQ A3 (GET A3 'ARRAY) B (GET B 'ARRAY) X (GET X 'ARRAY))
  580. ;;     (setq  AA (*ARRAY nil 'T (f1+ N) (f1+ N)))
  581. ;;     (setq  ROW (*ARRAY nil 'T (f1+ N)))
  582. ;;     (SETQ FLAGINV (= M 0.))
  583. ;;     (COND (FLAGINV (SETQ M N)))
  584. ;;     (DO ((I 1. (f1+ I)))
  585. ;;         ((> I N))
  586. ;;         (DO ((J 1. (f1+ J)))
  587. ;;         ((> J N))
  588. ;;         (STORE (AA I J)
  589. ;;                (NOT (EQUAL (AREF A3 I J) '(0. . 1.))))))
  590. ;;     (DO ((I 1. (f1+ I)))
  591. ;;         ((> I N))
  592. ;;         (DO ((J 1. (f1+ J)))
  593. ;;         ((> J M))
  594. ;;         (STORE (AREF X I J)
  595. ;;                (COND (FLAGINV (EQ I J))
  596. ;;                  (T (EQUAL (AREF B I J)
  597. ;;                        '(0. . 1.)))))))
  598. ;;     (DO ((I 1. (f1+ I))) ((> I N)) (STORE (ROW I) I))
  599. ;;         ;FORWARD ELIMINATION.
  600. ;;     (DO ((I 1. (f1+ I)))
  601. ;;         ((> I N))
  602. ;;         (SETQ FLAG-NONSINGULAR
  603. ;;           (DO ((II I (f1+ II)))
  604. ;;               ((> II N) NIL)
  605. ;;               (COND ((AA (ROW II) I)
  606. ;;                  (TMEXCHANGE 'ROW II I)
  607. ;;                  (RETURN T)))))
  608. ;;         (COND ((NULL FLAG-NONSINGULAR) (RETURN NIL)))
  609. ;;         (DO ((II (f1+ I) (f1+ II)))
  610. ;;         ((> II N))
  611. ;;         (COND ((AA (ROW II) I)
  612. ;;                (DO ((JJ (f1+ I) (f1+ JJ)))
  613. ;;                ((> JJ N))
  614. ;;                (STORE (AA (ROW II) JJ)
  615. ;;                   (OR (AA (ROW I) JJ)
  616. ;;                       (AA (ROW II) JJ))))
  617. ;;                (DO ((JJ 1. (f1+ JJ)))
  618. ;;                ((> JJ M))
  619. ;;                (STORE (AREF X (ROW II) JJ)
  620. ;;                   (OR (AREF X (ROW I) JJ)
  621. ;;                       (AREF X (ROW II) JJ))))))))
  622. ;;     (COND ((NULL FLAG-NONSINGULAR) (GO EXIT)))       ;GET OUT  BACKWARD SUBSTITUTION
  623. ;;     (DO ((I (f1- N) (f1- I)))
  624. ;;         ((< I 1.))
  625. ;;         (DO ((L 1. (f1+ L)))
  626. ;;         ((> L M))
  627. ;;         (STORE (AREF X (ROW I) L)
  628. ;;                (OR (AREF X (ROW I) L)
  629. ;;                (DO ((J (f1+ I) (f1+ J)) (SUM))
  630. ;;                    ((> J N) SUM)
  631. ;;                    (SETQ SUM
  632. ;;                      (OR SUM
  633. ;;                      (AND (AA (ROW I) J)
  634. ;;                           (AREF
  635. ;;                              X
  636. ;;                              (ROW J)
  637. ;;                              L)))))))))
  638. ;;            ;RECOVER THE ORDER.
  639. ;;     (TMPERMUTE 'X N M 0. 0. 'ROW N 'ROW)
  640. ;;    EXIT (*TMREARRAY 'ROW) (*TMREARRAY 'AA) (RETURN FLAG-NONSINGULAR)))
  641.  
  642. ;TMPERMUTE PERMUTES THE ROWS OR COLUMNS
  643. ;OF THE N*M MATRIX AX ACCORDING TO THE
  644. ;SPECIFICATION OF INDEXLIST. THE FLAG
  645. ;MUST BE SET 'ROW IF ROW PERMUTATION IS
  646. ;DESIRED , OR 'COL OTHERWISE. THE RESULT
  647. ;IS IN AX. NM IS THE DIMENSION OF
  648. ;INDEXLIST.
  649.  
  650. (DEFUN TMPERMUTE (AX N M RBIAS CBIAS INDEXLIST NM FLAG) 
  651.        (PROG (K L) 
  652. ;         (SETQ AX (GET AX 'array) 
  653. ;           INDEXLIST (GET INDEXLIST 'array))
  654.          (setq ax (get-array-pointer ax))
  655.          (setq indexlist (get-array-pointer indexlist))
  656.          (ARRAY *INDX* T (f1+ NM))
  657.          (DO ((I 1. (f1+ I)))
  658.          ((> I NM))
  659.          (STORE (aref *INDX* I) (AREF INDEXLIST I)))
  660.          (DO ((I 1. (f1+ I)))
  661.          ((> I NM))
  662.          (COND ((NOT (= (aref *INDX* I) I))
  663.             (PROG NIL 
  664.                   (TMMOVE AX N M RBIAS CBIAS I 0. FLAG)
  665.                   (SETQ L I)
  666.              LOOP (SETQ K (aref *INDX* L))
  667.                   (STORE (aref *INDX* L) L)
  668.                   (COND ((= K I)
  669.                      (TMMOVE AX
  670.                          N
  671.                          M
  672.                          RBIAS
  673.                          CBIAS
  674.                          0.
  675.                          L
  676.                          FLAG))
  677.                     (T (TMMOVE AX
  678.                            N
  679.                            M
  680.                            RBIAS
  681.                            CBIAS
  682.                            K
  683.                            L
  684.                            FLAG)
  685.                        (SETQ L K)
  686.                        (GO LOOP)))))))
  687.          (*TMREARRAY '*INDX*))) 
  688.  
  689. (DEFUN TMMOVE (AX N M RBIAS CBIAS I J FLAG) 
  690.        (PROG (LL)
  691.          (setq ax (get-array-pointer ax))
  692.          (SETQ LL (COND ((EQ FLAG '*ROW*) (f- M CBIAS))
  693.                 (T (f- N RBIAS))))
  694.          (DO ((K 1. (f1+ K)))
  695.          ((> K LL))
  696.          (COND ((EQ FLAG '*ROW*)
  697.             (STORE (AREF
  698.                       AX
  699.                       (f+ RBIAS J)
  700.                       (f+ CBIAS K))
  701.                    (AREF
  702.                       AX
  703.                       (f+ RBIAS I)
  704.                       (f+ CBIAS K))))
  705.                (T (STORE (AREF
  706.                         AX
  707.                         (f+ RBIAS K)
  708.                         (f+ CBIAS J))
  709.                  (AREF
  710.                         AX
  711.                         (f+ RBIAS K)
  712.                         (f+ CBIAS I))))))))
  713.  
  714. ;TMSYMETRICP CHECKS THE SYMETRY OF THE MATRIX.
  715.  
  716. (DEFUN TMSYMETRICP        (A3 N)
  717.        (SETQ A3 (GET-array-pointer A3))
  718.        (DO ((I 1. (f1+ I)))
  719.        ((> I N) T)
  720.        (COND ((NULL (DO ((J (f1+ I) (f1+ J)))
  721.                 ((> J N) T)
  722.                 (COND ((NOT (EQUAL (AREF
  723.                               A3
  724.                               I
  725.                               J)
  726.                            (AREF
  727.                               A3
  728.                               J
  729.                               I)))
  730.                    (RETURN NIL)))))
  731.           (RETURN NIL)))))
  732.  
  733. ;TMLATTICE CHECKS THE "LATTICE"
  734. ;STRUCTURE OF THE MATRIX A. IT RETURNS
  735. ;NIL IF THE MATRIX IS "OBVIOUSLY"
  736. ;SINGULAR. OTHERWISE IT RETURNS A LIST
  737. ;(L1 L2 ... LM) WHERE M IS THE NUMBER OF
  738. ;BLOCKS (STRONGLY CONNECTED SUBGRAPHS),
  739. ;AND L1 L2 ... ARE LIST OF ROW AND
  740. ;COLUMN NUBERS WHICH BELONG TO EACH
  741. ;BLOCKS. THE LIST LOOKS LIKE ((R1 C1)
  742. ;(R2 C2) ...) WHERE R R'S ARE ROWS AND
  743. ;C'S ARE COLUMMS.
  744.  
  745. (DEFUN TMLATTICE (A3 XROW XCOL N) 
  746.        (PROG (RES) 
  747.          (setq a3 (get-array-pointer a3))
  748.          (setq xrow (get-array-pointer xrow))
  749.          (setq xcol (get-array-pointer xcol))
  750.          (setq *b* (*ARRAY nil T (f1+ N) (f1+ N)))
  751.          (setq *ROW* (*ARRAY nil T (f1+ N)))
  752.          (setq *col* (*ARRAY nil  T (f1+ N)))
  753.          (DO ((I 1. (f1+ I)))
  754.          ((> I N))
  755.          (DO ((J 1. (f1+ J)))
  756.              ((> J N))
  757.              (STORE (aref *B* I J)
  758.                 (NOT (EQUAL (AREF A3 I J)
  759.                     '(0. . 1.))))))
  760.          (DO ((I 0. (f1+ I)))
  761.          ((> I N))
  762.          (STORE (aref *ROW* I) I)
  763.          (STORE (aref *COL* I) I))
  764.          (COND ((NULL (TMPIVOT-ISOLATE 1.))
  765.             (SETQ RES NIL)
  766.             (GO EXIT)))
  767.          (DO ((I 1. (f1+ I)))
  768.          ((> I N))
  769.          (STORE (aref *B* (aref *ROW* I) (aref *COL* I)) I)
  770.          (STORE (aref *B* (aref *ROW* I) (aref *COL* 0.)) T
  771.             ))
  772.          (TMLATTICE1 1.)
  773.          (SETQ RES (TMSORT-LATTICE XROW XCOL))
  774.     EXIT (*TMREARRAY '*B*)
  775.          (*TMREARRAY '*ROW*)
  776.          (*TMREARRAY '*COL*)
  777.          (RETURN RES))) 
  778.  
  779. (DEFUN TMLATTICE1 (K) 
  780.        (COND ((= K N) NIL)
  781.          (T (TMLATTICE1 (f1+ K))
  782.         (DO ((LOOPPATH))
  783.             (NIL)
  784.             (COND ((SETQ LOOPPATH (TMPATHP K K))
  785.                (TMUNIFY-LOOP K (CDR LOOPPATH)))
  786.               (T (RETURN NIL))))))) 
  787.  
  788. (DEFUN TMPATHP (J K) 
  789.        (COND ((EQUAL (aref *B* (aref *ROW* J) (aref *COL* K)) T) (LIST J K))
  790.          (T (DO ((JJ K (f1+ JJ)) (PATH))
  791.             ((> JJ N))
  792.             (COND ((AND (EQUAL (aref *B* (aref *ROW* J) (aref *COL* JJ)) T)
  793.                 (SETQ PATH (TMPATHP JJ K)))
  794.                (RETURN (CONS J PATH)))))))) 
  795.  
  796. (DEFUN TMUNIFY-LOOP (K CHAIN) 
  797.        (PROG (L DUMMYK DUMMYL) 
  798.          (SETQ L (CAR CHAIN))
  799.          (COND ((= L K) (RETURN NIL)))
  800.          (SETQ DUMMYK (aref *B* (aref *ROW* K) (aref *COL* K)))
  801.          (SETQ DUMMYL (aref *B* (aref *ROW* L) (aref *COL* L)))
  802.          (STORE (aref *B* (aref *ROW* K) (aref *COL* K)) NIL)
  803.          (STORE (aref *B* (aref *ROW* L) (aref *COL* L)) NIL)
  804.          (DO ((I 1. (f1+ I)))
  805.          ((> I N))
  806.          (STORE (aref *B* (aref *ROW* K) (aref *COL* I))
  807.             (OR (aref *B* (aref *ROW* K) (aref *COL* I)) (aref *B* (aref *ROW* L) (aref *COL* I))))
  808.          (STORE (aref *B* (aref *ROW* I) (aref *COL* K))
  809.             (OR (aref *B* (aref *ROW* I) (aref *COL* K)) (aref *B* (aref *ROW* I) (aref *COL* L))))
  810.          (STORE (aref *B* (aref *ROW* L) (aref *COL* I)) NIL)
  811.          (STORE (aref *B* (aref *ROW* I) (aref *COL* L)) NIL))
  812.          (STORE (aref *B* (aref *ROW* K) (aref *COL* K)) DUMMYL)
  813.          (STORE (aref *B* (aref *ROW* L) (aref *COL* L)) DUMMYK)
  814.          (STORE (aref *B* (aref *ROW* K) (aref *COL* 0.)) T)
  815.          (STORE (aref *B* (aref *ROW* L) (aref *COL* 0.)) NIL)
  816.          (TMUNIFY-LOOP K (CDR CHAIN)))) 
  817.  
  818. (DEFUN TMSORT-LATTICE (XROW XCOL) 
  819.        (PROG (NODELIST RESULT) 
  820.          (SETQ NODELIST (TMSORT1))
  821.          (SETQ 
  822.           RESULT
  823.           (DO ((X NODELIST (CDR X)) (RESULT))
  824.           ((NULL X) RESULT)
  825.           (SETQ RESULT
  826.             (CONS (DO ((NEXT (aref *B* (aref *ROW* (CAR X))
  827.                         (aref *COL* (CAR X)))
  828.                      (aref *B* (aref *ROW* NEXT) (aref *COL* NEXT)))
  829.                    (RES))
  830.                   ((= NEXT (CAR X))
  831.                    (CONS (LIST (aref *ROW* NEXT) (aref *COL* NEXT))
  832.                      RES))
  833.                   (SETQ RES
  834.                     (CONS (LIST (aref *ROW* NEXT)
  835.                             (aref *COL* NEXT))
  836.                           RES)))
  837.                   RESULT))))
  838.          (DO ((LIST1 RESULT (CDR LIST1)) (I 1.))
  839.          ((NULL LIST1))
  840.          (DO ((LIST2 (CAR LIST1) (CDR LIST2)))
  841.              ((NULL LIST2))
  842.              (STORE (AREF XROW I) (CAAR LIST2))
  843.              (STORE (AREF XCOL I) (CADAR LIST2))
  844.              (SETQ I (f1+ I))))
  845.          (RETURN RESULT))) 
  846.  
  847. ;; (DEFUN TMLESS (I J) (B (ROW I) (COL J))) 
  848.  
  849. (DEFUN TMSORT1 NIL 
  850.        (DO ((I 1. (f1+ I)) (RESULT))
  851.        ((> I N) RESULT)
  852.        (COND ((AND (aref *B* (aref *ROW* I) (aref *COL* 0.)) (TMMAXP I))
  853.           (DO ((J 1. (f1+ J)))
  854.               ((> J N))
  855.               (COND ((NOT (= J I))
  856.                  (STORE (aref *B* (aref *ROW* I) (aref *COL* J)) NIL))))
  857.           (STORE (aref *B* (aref *ROW* I) (aref *COL* 0.)) NIL)
  858.           (SETQ RESULT (CONS I RESULT))
  859.           (SETQ I 0.))))) 
  860.  
  861. (DEFUN TMMAXP (I) 
  862.        (DO ((J 1. (f1+ J)))
  863.        ((> J N) T)
  864.        (COND ((AND (NOT (= I J)) (aref *B* (aref *ROW* J) (aref *COL* I)))
  865.           (RETURN NIL)))))
  866.  
  867. ;;UNPIVOT IS USED IN PAUL WANG'S PROGRAM
  868. ;;TO RECOVER THE PIVOTTING. TO GET THE
  869. ;;INVERSE OF A, PAUL'S PROGRAM COMPUTES
  870. ;;THE INVERSE OF U*A*V BECAUSE OF
  871. ;;BLOCKING. LET THE INVERSE Y. THEN
  872. ;;A^^-1=V*Y*U. WHERE U AND V ARE
  873. ;;FUNDAMENTAL TRANSFORMATION
  874.  ;;(PERMUTATION). UNPIVOT DOES THIS,
  875. ;;NAMELY, GIVEN A MATRIX A3, INDEX ROW
  876. ;;AND COL ,WHICH CORRESPONDS TO THE Y , U
  877. ;; AND V, RESPECTIVELY, IT COMPUTES V*Y*U
  878. ;;AND RETURNS IT TO THE SAME ARGUMENT A.
  879.  
  880. (DEFUN TMUNPIVOT (A3 *ROW* *COL* N M) 
  881.        (PROG NIL 
  882.          (setq *col* (get-array-pointer *col*))
  883.          (setq *row* (get-array-pointer *row*))
  884.          (setq *rowinv* (*ARRAY nil T (f1+ N)))
  885.          (setq *colinv* (*ARRAY nil T (f1+ N)))
  886.          (DO ((I 1. (f1+ I)))
  887.          ((> I N))
  888.          (STORE (aref *ROWINV* (AREF *ROW* I)) I))
  889.          (DO ((I 1. (f1+ I)))
  890.          ((> I N))
  891.          (STORE (aref *COLINV* (AREF *COL* I)) I))
  892.          (TMPERMUTE A3 N M 0. N '*COLINV* N '*ROW*)
  893.          (TMPERMUTE A3 N M 0. N '*ROWINV* N '*COL*)
  894.          (*TMREARRAY '*ROWINV*)
  895.          (*TMREARRAY '*COLINV*))) 
  896.  
  897. #-NIL
  898. (DECLARE-TOP(UNSPECIAL N  #-cl vlist NX IX))
  899.